TestMAPlot

library(Glimma)
library(limma)
library(GlimmaV2)

data(lymphomaRNAseq)
rnaseq <- lymphomaRNAseq

# 7 samples
# each instance is a gene
# attribute is number of counts per sample
# data.frame(rnaseq$counts)

# add lane
groups <- data.frame(genotype=rnaseq$samples$group,
                     lane= as.character(c(rep(4,5),3,3)),
                     miscCont=c(rep(4000,5),300,250),
                     miscDisc=c("blue","red",rep("green",5)))

# add libsize
groups <- rnaseq$samples$group

# fit
design <- model.matrix(~0+groups)
contrasts <- cbind(Smchd1null.vs.WT=c(-1,1))

# convert raw counts to logCPM values by automatically extracting libsizes and normalisation factors from x
vm <- voomWithQualityWeights(rnaseq, design=design)
fit <- lmFit(vm, design=design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
dt <- decideTests(fit)
summary(dt)

##        Smchd1null.vs.WT
## Down               1478
## NotSig            24424
## Up                 1277

testing <- function(fit, counts, groups=1:ncol(counts), logTransform=FALSE)
{
  
  if (is.null(colnames(counts)))
  {
    samples <- seq_len(ncol(counts))
  }
  else
  {
    samples <- colnames(counts)
  }
  
  # process counts
  cols <- colnames(counts)
  rownames(counts) <- make.names(cols)
  if (logTransform) 
  {
      counts <- as.matrix(edgeR::cpm(counts, log=TRUE))
  } 
  else 
  {
      counts <- as.matrix(counts)
  }
  counts <- t(counts)
  return(counts)
  if (!is.numeric(groups)) 
  {
      # reorder samples to group levels
      groups <- factor(groups)
      counts <- counts[order(groups), ]
      sample.cols <- sample.cols[order(groups)]
      samples <- samples[order(groups)]
      groups <- sort(groups)
  }

  expData <- data.frame(
      Sample = samples,
      cols = as.hexcol(sample.cols),
      Group = groups,
      counts)
  
  expData

}

# object has a bunch of genes
nrow(fit$genes)

## [1] 27179

# fold change is the attribute here
# each instance is a gene;
# the col vector names are the geneIDs
# the col vector attribute is fold change (Smchd1null vs WT)
# data.frame(fit$coefficients)

# log cpm is the attribute here
# averaged across all arrays in the original linear model fit
# each instance is a gene;
# the col vector names are the geneIDs
# the col vector attribute is mean expression
# data.frame(fit$Amean)


# for now assume fit object

# this is where info for the MD scatter plot comes from!
# data.frame(vm$E)
# fit$design

GlimmaV2("MA", fit)

## Warning in readLines(con): incomplete final line found on 'C:/Users/hkari/
## OneDrive/onenote/R/win-library/3.6/GlimmaV2/htmlwidgets/GlimmaV2.yaml'